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Ab initio multi-reference configuration interaction calcu- 
lations are performed for the Si(100) surface using a cluster 
approach. The convergence with respect to the cluster size is 
checked and the final results are taken from a Si32H2g cluster 
which models two dimers and six bulk layers. We find for the 
ideal as well as for the p(l x 2) reconstruction a singlet ground 
state consisting of several configurations. The energy gain due 
to forming the symmetric dimer in the p(l x 2) structure is 
1.75 eV, the bond length of the dimer is 2.35 A which is very 
close to the bulk value. In contradiction to the LDA results 
and in agreement with previous correlation calculations we do 
not find an asymmetric p(l x 2) structure. 
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I. INTRODUCTION 

The Si(100) surface, is one of the bes!— examined, 
both experimentallyErLl and theoreticallyEro, semicon- 
ductor surfaces due to its great technological imparlance 
for semiconductor devices. JJEED experiments™!, X- 
ray diffraction nieasurementsu'Ll, scanning tunneling mi- 
croscopy (STM)H and surface core-level photoelectron 
diffraction experimental find down to temperatures of 
90K that on the Si(100) surface alternating buckled 
dimers in p(2x2) or c(4x2) structure are formed. These 
experimental findings are in agreement with those of var- 
ious density functional calculations (DFT) applying the 
local deiis^y-apjproximation (LDA) to a slab model of the 
surfaceU'Ej'Oia. Applying the LDA to cluster models of 
this surfaceEJcH the energy difference between the buck- 
led reconstruction and a reconstruction forming symmet- 
ric dimers is very small, so they can not draw a conclud- 
ing result on the ground state structure. Earlier work 
done with semi-empirical quantum chemical methodstl 
favours the asymmetricjeconstruction as a closed-shell 
Hartree-Fock treatments does. Including electronic cor- 
relations at the generalized valence bond leveB or within 
a configuration interaction (CI) approachEJ yields a sym- 
metrically dimerized ground state similar-to the findings 
of unrestricted Hartree-Fock calculationsEl 
The wide spread density functional methods rely on the 
ground-state density and do not provide a many-body 
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FIG. 1. The figures show sketches of the orbital structure 
of the ideal surface (a) , and of a closed-shell (c) and open shell 
(b) structure of the reconstructed surface. 



wave function. The electronic correlations are included 
implicitly only via the exchange-correlation functional. 
Regarding the process forming a dimer on the surface the 
question arises whether DFT can cope with bond break- 
ing and bond formation. In the ideal Si(100) surface 
each surface atom has two singly occupied sp 3 hybrids 
(see Fig.^ a). In a simple-minded picture the surface 
energy is lowered by reducing the number of singly oc- 
cupied orbitals. A bond is formed between two surface 
atoms yielding the (1x2) dimer structure. The nature of 
the bond can have two limiting cases: If a double-bond 
is generated between the two Si atoms (see Fig.^Jb), this 
would lead to a closed-shell structure of the surface. Or, 
a single-bond is generated leaving on each surface atom a 
singly occupied orbital (see Fig.^jc). The true electronic 
structure of the surface can only be solved within a multi- 
determinant approach, which can describe all configura- 
tions between the two limiting cases. A second advantage 
of this approach is that the ideal and reconstructed sur- 
face are described at the same level of accuracy because 
we apply no restrictions on the wave function as single- 
reference approaches do. Up to now the multi-reference 
methods were applied to adsorption and diffusion pro- 
cesses of hydrogen on Si(100)t3. The most sophisticated 
correlation method apnlLsd to the reconstruction itself 
is by Jing and WhittentL-J who made a two-determinant 
ansatz at the MCSCF level and perform a CI calculation 
using this reference. 

We want to present a multi-reference self-consistent field 
(MCSCF) calculation for the ideal and reconstructed 
Si(100) surface to cope with the difficult bond structure. 
To deal with the dynamical correlations we supplement it 
with a multi-reference configuration interaction (MRCI) 
calculation. In addition we perform LDA and gradient 
corrected (GGA) density functional calculations for a di- 
rect comparison. 

The paper is organized as follows: In Sec. II we present 
the multi-reference methods and provide the technical de- 
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tails. The resulting wave function is discussed in Sec. III. 
In Sec. IV we present the reconstruction energy, the ge- 
ometry of the symmetric dimer and the question of dimer 
buckling. Conclusion follows in Sec.V. 

II. MULTI-REFERENCE METHODS 

Only closed-shell systems or high-spin open-shell sys- 
tem can be described within a single-determinant ap- 
proach such as Hartree-Fock or single-reference con- 
figuration interaction (CI) methods to include correla- 
tions. In all other cases the suitable method is a multi- 
configuration self-consistent field (MCSCF) approach. 
The wave function is written as 

I^mcscf >= X) C 'I* J >' W 

where > is a determinant built up from molecular or- 
bitals and cj are the so-called Cl-coefficients. The differ- 
ence to a single- reference CI approach is that not only the 
Cl-coefficients are optimized but also the coefficients of 
the orthonormal orbitals building up the >. There- 
fore the number of determinants involved in the MCSCF 
should be small compared with conventional CI. In the 
case of the Si(100) surface we select all possible deter- 
minants combining the four singly occupied sp 3 hybrids 
to a singlet. In addition we reoptimize the closed-shell 
orbitals of the first layer to account for the influence of 
the reconstruction on the bulk. The electronic correla- 
tions treated by a MCSCF ansatz are often called static 
correlations. 

To describe the dynamical correlations we apply on top 
of the MCSCF a multi-reference configuration interaction 
(MRCI) with single and double excitations. 

I*MRCI >= 1 + XX C i C " + J2 C C H C M C * l*MCSCF 

(2) 

where fj, and v numbers the occupied orbitals and i and 
j the unoccupied ones. The coefficients 7jf and rj^ are 
determined variationally in the MRCI procedure. We 
restrict the dynamically correlated space to the surface 
orbitals. 

We perform all calculations with, the quantum-chemical 
program package MOLPR094B, only for thft_pFT- 
calculations we use the latest version M0LPR096O. The 
smallest possible cluster with one dimer is Sig, the bulk- 
like bonds are saturated with 12 hydrogens. To test the 
influence of the bulk layers we perform calculations on 
Sii7H2o and Si3iH32, where the number of bulk layers 
increases. The Si — Si distance in the bulk is taken from 
experiment (2.352A), the distance to the saturating hy- 
drogen atoms is chosen equal to the one in SiH4 (1.480A). 




FIG. 2. The figure shows all Si^Hy clusters which are used 
in the calculations; the surface atoms are shaded; the hydro- 
gen atoms are not drawn. 

In test calculations the geometry of the symmetric re- 
construction is chosen such that all bonds have the bond 
length of the bulk. To test the interaction of the dimers 
we consider a Sii5H 16 cluster which contains two par- 
allel dimers and Si23H24 which has two dimers in line. 
Additionally we used a Si2iH 2 o cluster with three paral- 
lel dimers. The final results are obtained from Si32H28, 
which simulates six bulk layers and two dimers. All clus- 
ters are shown in Fig. |^. 

We apply a 4-valence-electron scalar-relativistic energy- 
consistent pseudopoteotial for SiE3 with the correspond- 
ing valence basis setEE It is supplemented by various 
sets of polarization functions, which are optimized by 
performing correlated calculations for the free atom. For 
hydrogen we apply a minimal basis consisting of the fully 
contracted Is functions of DunningH. 

III. MANY-BODY WAVE-FUNCTION OF THE 
IDEAL AND RECONSTRUCTED SURFACE 

For the ideal surface we perform a high-spin (S = 2) 
single-determinant Hartree-Fock calculation to generate 
the starting orbitals for the MCSCF calculation in which 
the four singly occupied orbitals in the active space are 
combined to a singlet ground-state. In addition the 
closed-shell orbitals of the first bulk layer are reopti- 
mized. One way to discuss the resulting orbitals is to 
calculate the natural orbitals which diagonalize the den- 
sity matrix (see Table ||) . We find three quite close-lying 
orbitals two of which have binding and one has antibind- 
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TABLE I. The mean occupation number of the natural 



orbitals and their binding character is shown both for the 
ideal and reconstructed surface. 
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71 


ideal 


1. bind. 


1.52 




2. bind. 


1.06 




3. antib. 


0.94 




4. antib. 


0.48 


sym. reconstr. 


1. bind. 


1.98 




2. bind. 


1.63 




3. antib. 


0.37 




4. antib. 


0.02 



ing character. Even the second antibinding one has a 
mean occupation of half an electron. Regarding the co- 
efficients of the MCSCF expansion there are ten config- 
urations with a coefficient greater than 0.1. This shows 
clearly that the ideal surface can not be described within 
a single-determinant approach. The singlet closed-shell 
HF energy lies about 1 eV higher in energy than the cor- 
responding MCSCF value. 

For the symmetric reconstructed surface we start the 
MCSCF calculation with the orbitals of a triplet HF 
determinant. But the active space of the MCSCF cal- 
culation is exactly the same as in the ideal case. The 
natural orbitals look quite different. A low- lying binding 
orbital is doubly occupied forming a single bond between 
the surface atoms. A second binding orbital is occupied 
with 1.6 electrons. A double occupancy of this orbital 
would correspond to a double-bond between the surface 
atoms. The significant population of the first antibind- 
ing orbital with 0.4 electrons contradicts the closed-shell 
picture. The MCSCF expansion is dominated by the con- 
figuration where the two binding orbitals are occupied, 
but the occupation of the first antibinding orbital is not 
negligible. In a single-determinant description the triplet 
state would be energetically favoured by 0.9 eV leading 
to an unphysical ferromagnetic surface. In the MCSCF 
calculation the triplet state is about 0.3 eV higher in en- 
ergy than the singlet ground state. 

In both cases, i.e., in the ideal and in the symmetrically 
reconstructed surface we find a singlet ground state which 
can be only described within a multi-configuration ap- 
proach. Through this we can provide for both cases an 
adequate description. So we can compare the ground- 
state energies of the two structures yielding the recon- 
struction energy. 

IV. RESULTS AND DISCUSSION 

A. Basis set dependence 

In a first step we want to test the quality of the ba- 
sis set used. We calculate the reconstruction energy on 
the MCSCF and MRCI level for the Si 9 Hi 2 cluster. Al- 



TABLE II. Influence of the atomic basis set on the recon- 
struction energy at different correlation levels. The calcula- 
tions are performed in the Si9Hi2 cluster. 



basis 


ASmcscf [eV] 


A£ M RCi[eV] 


[2s2p] 


0.90 


1.08 


red. basis 


1.17 


1.54 


[2s2p]ld 


1.26 


1.56 


(4s4p)ld 


1.27 


1.59 


[3s3p]2d 


1.29 


1.60 


[3s3p]2dl/ 


1.33 


1.73 


(4s4p)3d2/lg 


1.34 


1.73 



though the question of the basis is a technical one, it is 
possible to draw some physical conclusions. 
The default basis is of valence-double-C quality 
(AsAp) I [2s2p] supplemented by one d polarization func- 
tion with exponent 0.4. In Table || the results are listed. 
Calculating the reconstruction energy without any d- 
function yields only 70% of the value with the default 
basis. Applying a d-function only for the surface and 
the first layer atoms (reduced basis) we determine 93% 
of the default MCSCF energy and 97% of the default 
MRCI energy. The better value of the MRCI is clear, 
because in this calculation only the surface orbitals are 
involved, whereas in the MCSCF the first layer is re- 
optimized, too. This basis set will be our choice for the 
larger clusters. The increase of the reconstruction energy 
when supplementing a d polarisation function indicates 
the important role of the electronic correlations on the 
formation of the bond in the dimer. Uncontracting to- 
tally the valence basis set yields only a small increase 
in the reconstruction energy. This shows that the [2s2p] 
basis is of equal quality for both geometries. The same 
holds true if we split the d-exponent into two (0.23, 0.8). 
To receive a balanced basis we uncontract the valence 
basis to [3s3p]. Applying an additional /-function (ex- 
ponent 0.35) the reconstruction energy rises by 6% on 
the MCSCF level and by 10% on the MRCI level. This 
indicates that to correlate the 7r-bond the /-function is 
more important in the reconstructed geometry, especially 
for the dynamical correlation covered by the MRCI. Un- 
contracting the valence basis set totally to (4s4p) and 
adding a (3d2flg) polarization set (exponents from Ref.[ 
p4| ) yields nearly no changes (smaller than 1%) on the 
reconstruction energy on the MCSCF and MRCI level. 
This indicates that the basis set limit is reached for the 
reconstruction energy . 

B. Influence of bulk layers 

To check the influence of the bulk layers on the sur- 
face reconstruction we perform MCSCF and MRCI cal- 
culations for three different clusters using the reduced 
basis. SigHi2 simulates three bulk layers, Sii7H2o four 
layers and Si3iH32 five layers whereas all have one sur- 
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TABLE III. Influence of the bulk layers on the reconstruc- 
tion energy at different correlation levels. The calculations 
are performed in the reduced basis. 



cluster 


ABmcscf [eV] 


A£ M RCi[eV] 


SigHi2 


1.17 


1.54 


Sii7H2o 


1.33 


1.67 


Si3lH32 


1.30 


1.61 



TABLE IV. Influence of the dimer interaction on the re- 
construction energy at different correlation levels. The calcu- 
lations are performed in the reduced basis. 



cluster 


A£ M cscF[eV] 


AS MR ci[eV] 


A£ LDA [eV] 


AS GGA [eV 


SigHi2 


1.17 


1.54 


2.01 


1.43 


Sii5Hi 6 


1.15 


1.46 


1.85 


1.27 


Si2iH2o 


1.13 




1.76 


1.27 


Si23H24 


1.15 


1.46 







face dimer; for results see Table III. The changes between 
three and four layers are quite significant, i.e. up to 14% 
of the reconstruction energy at the MCSCF level. This 
shows that the SigHi 2 cluster is not fully reliable for de- 
scribing the surface. The difference between four and five 
layers is up to 4% and is tolerable for the determination 
of the reconstruction energy. But in our opinion it is not 
possible to calculate reliable relaxed distances of the bulk 
atoms, which is a shortcoming of the cluster approach. 
Experimentally it is provcrfl that there are relaxations 
up to the fifth bulk layer. 



C. Interaction of dimers 

Although we have checked in the last section, that the 
3-layer clusters can not fully describe the surface recon- 
struction we select four different 3-layer clusters in order 
to test the interaction of the dimers. SigHi2 has one 
dimer, Sii5rl 16 has two parallel dimers and Si2iH 2 o has 
three. Si23H 2 4 has two dimers in line. Concerning the 
reconstruction energy on the MCSCF level, all selected 
clusters differ by at most 4% concluding that the interac- 
tion on the MCSCF level is weak (see Table 0). On the 
MRCI level we calculate only the cluster with one and 
two parallel dimers. The difference is larger than at the 
MCSCF level indicating that the dynamical correlations 
are more important for the dimer interaction than are 
the static correlations. 

Treating the clusters in a LDA approach (Slater- 
Dirac exchangeo and Vosko-Wilk-Nusair correlation 
functionaEj) the interaction between the dimers is larger 
by up to 12%. This can be understood from the wave 
function structure resulting from the LDA. If more than 
one dimer exist, there is an energy gain due to derealiza- 
tion of the surface electrons. In a gradient corrected den- 
sity functional approach (GGA, Beckje exchangee!! and 
Lee- Yang-Parr correlation functionalEj) we found a sim- 



ilar behaviour, but there the effect of derealization is 
smaller than in the LDA and occurs only between two 
dimers. This derealization energy is larger on the ideal 
surface than on the reconstructed surface due to the close 
lying and therefore well mixing states. This results in an 
overestimation of the dimer interaction on the LDA level. 
If the surface is treated on the MCSCF level the mixing 
between different states can occur within one dimer, so 
that both ideal and reconstructed surface are treated at 
the same level of accuracy. 



D. Surface reconstruction energy — 
different methods 



comparison of 



Combining the results of the previous sections we se- 
lect the Si3 2 H 2 § cluster for our final calculation. This 
cluster has two dimers and the bulk is modelled up to 
the sixth bulk layer. Due to the quite large size of the 
cluster we can treat it only with the reduced basis. The 
surface reconstruction energy at the MCSCF level is 1.13 
eV. Including dynamical correlations yields an increase 
of nearly 40% resulting in 1.56 eV. Concerning the lim- 
ited basis used in this calculation we can estimate based 
on our basis set test calculations a final value of 1.75 
eV. The error due to the cluster size is estimated to 
be by ±5% . Comparing our result-,with various LDA- 
slab, calculations (Dabrowskij-ejt al.Ej 1.5 eV; Fritsch et 
al.N 1.51 eV; Ramstad et al£3 1.8 eV; Roberts et al.til 
2.02 eV) the agreement is remarkably good. In order 
to compare the results of different methods in more de- 
tail we apply the different methods to the same cluster 
with the same basis set. We choose again the SigHi 2 
cluster using the [2s2p]ld basis. DFT methods are ap- 
plied with two different exchange correlation functions. 
The LDA described above yield a reconstruction energy 
of 2.01 eV which is nearly 30% above the MRCI value 
of 1.56 eV. This value is certainly lowered due to dimer 
interaction as described in the previous section. Using 
the GGA approach the reconstruction energy (1.43 eV) 
is underestimated by 8% . The DFT treatments yield 
a singlet ground state for the reconstructed surface in 
contradiction to a single-determinant Hartree-Fock cal- 
culation, but the reconstruction energy is dependent on 
the chosen functional. The results show the known fea- 
ture that LDA overestimates binding whereas the GGA 
gives a substantial improvement. 



E. Geometry of the symmetric reconstruction 

We argued in the previous sections that the clusters 
are too small for treating the relaxations of the bulk 
layers due to surface reconstruction. But it is reason- 
able to optimize the dimer bond length (6dimcr) and 
the bond length of the surface atoms to the first layer 
atoms (b\, layer); for results see Table [v| To check the 
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TABLE V. Optimized dimer bond length and first layer 
bond length of the symmetric reconstructed surface. 



method 


h rXi 

Odimcr \J\\ 


Ol. layer [AJ 


MCSCF 


Si 9 H 


12 


2.41 


2.36 


MRCI 


Si 9 H 


12 


2.37 


2.36 


MCSCF 


Sii7H2o 


2.40 


2.36 


MRCI 


Sii7H2o 


2.35 


2.36 


LDA 


[ 


11 




2.23 


2.29 




[ 


H 




2.23 


2.27 




[ 


16 




2.26 


2.31 




FIG. 3. The figure shows a projection of the SigHi2 cluster 
where the displacements concerned for the asymmetric buck- 
ling are indicated. 

cluster-size dependence we perform the calculations for 
the SigHi2 cluster as well as for the SiiyH^o cluster within 
the [2s2p]ld basis. The results show, that the bond 
lengths are only slightly cluster dependent. The changes 
are less than 0.02A. The dynamical correlations reduce 
the dimer bond length by about 2% . Compared with 
the dimer bond-length of the asymmetric dimer from dif=. 
ferent experimcn ts (Bullock et al.til 2.25A; Jedrecy et al.til 
2.32A: Aono et alE 2.4±0.ll; Holland et all 2.47A; Yang 
et alfl 2.54A; Felici et al.0 2.67A) we agree quite well 
with the second and third measurement. LDA-slab cal- 
culations (see Table [v|) yield a dimer bond length quite 
below the bulk value. 



F. Asymmetric reconstruction 

In this section we want to arise the question of the 
dimer buckling. Up to now we allowed only for the 
symmetric reconstruction which maintains the local C2 V - 
symmetry. We test the buckling by lowering the sym- 
metry in the SigHi2 cluster with the [2s2p]ld basis. We 
allow for a buckling angle O and a displacement of the 
dimer parallel to the surface (see Fig. ||). The dimer 
bond length is fixed to the value obtained in the op- 
timization of the symmetric case. In the MCSCF and 
MRCI treatment the symmetric dimer has the lowest en- 
ergy contrary to the DFT calculations where we find a 
buckling of 9° in LDA and of 8° in GGA. The energy 
gain is very small, i.e., 0.031 eV in LDA and 0.014 eV in 



GGA. From this we can conclude that the used approxi- 
mations to DFT overestimate the electronic contribution 
to the buckling due to an overestimation of the closed- 
shell structure. 

The experimentally observed buckling could be driven by 
relaxation of the bulk atoms. Although the SigHi2 cluster 
is quite small to resolve relaxations we performed some 
test calculations within this model. We allow for the 
first bulk layer atoms to relax in all directions, whereas 
we regard for the second bulk layer only a relaxation per- 
pendicular to the surface (see Fig. ||). On the MCSCF 
level, where all orbitals are reoptimized we do not find 
any buckling although there is a significant relaxation of 
the first bulk layer atoms in the dimer direction of about 
0.08A and of the second bulk layer atoms perpendicular 
to the surface of about 0.09A. The other relaxations are 
smaller by about one order of magnitude. From this we 
can conclude that in our cluster model the bulk relax- 
ations do not drive the asymmetric buckling. 
To account for the dimer interaction we regard the 
Sii5Hi 6 cluster with two dimers. We model the asym- 
metric p(lx2) structure as well as the p(2x2) structure 
with alternating buckled dimers. For fixed buckling angle 
the p(2x2) structure is lower in energy than the asym- 
metric p(lx2) structure, e.g., for an angle of 2° by 0.018 
eV per dimer and for an angle of 10° by 0.060 eV. But 
still the symmetric dimer is lowest in energy, also if we 
allow for the relaxations we determined in the SigHi2 
cluster. This shows that the interaction of two dimers 
is not enough to stabilize a buckled structure. This is 
in .contradiction to the LDA cluster results from Yang et 
al.E£l who find a stable p(2x2) structure in the SiisHig 
cluster. 

Our finding of a symmetric ground state is in agree- 
ment with all previous correlation calculations, which 
show that the proper treatment of the electronic correla- 
tions stabilize the symmetric reconstruction. We can ex- 
clude an asymmetric. p(lx2) structure, which LDA slab 
calculationsEllB-Blia find about O.leV more stable than 
the symmetric one. As we can only calculate the interac- 
tion of two dimers, there is the possibility that the inter- 
action of many dimers can yield a buckled grouncLstate. 
But regarding the STM experiment of Badt et al.cl, who 
find a buckled c(2x4) structure at 90K, they claim that 
the buckling can be induced by a defect or a step at the 
surface. The defect-free surface is difficult to prepare, so 
the question of buckling of the Si(100) surface is experi- 
mentally still not solved. 

To address the question of metallicity of the symmetric 
reconstructed Si(100) surface, which is found in various 
DFT calculations and discussed in detail in Ref.[ p9[ , in 
a many-body description as applied here there can ex- 
ist a gap due to electronic correlations in the symmet- 
ric reconstruction, too. But in principle the question of 
metallicity cannot be solved in a cluster approach. 
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V. CONCLUSION 

We have performed MCSCF calculations for the ideal 
and reconstructed Si(100) surface. In both cases the 
many-body ground state is a singlet, but it can only 
be described by several determinants. On the recon- 
structed surface a dimcr is formed. The bond structure 
is a mixture between a double-bond and a single-bond 
with two singly occupied orbitals. To account for the 
dynamical correlations we apply on top of the MCSCF 
result a MRCI calculation. Determining the reconstruc- 
tion energy we find good agreement with various LDA 
slab calculations. The dimcr bond length agrees better 
with experiment than the LDA one. The asymmetric re- 
construction is overestimated in the LDA while our fully 
correlated treatment finds the symmetric dimcr more sta- 
ble regarding the electronic structure of one dimer and 
relaxations of the bulk atoms. The buckling observed in 
experiment can be due to a small concentration of defects 
on the Si(100) surface. It is possible that an alternating 
buckled structure is induced by defects and persist down 
to low temperatures due to dimer interactions. 
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